function varargout = cgparse0(varargin)

% CGPARSE0 Solve for PARSE parameters using a variant of CG
%
% CGPARSE0(FID,CGP,MODEL)
%
% Pass data parameters in with the structure FID, CG algorithm control
% parameters with the structure CGP, and initial model parameter estimates
% with MODEL.  Output is returned in a model structure.  If coded correctly,
% the model should be determined from the initial estimate structure (future work).
%
% CGPARSE0(FID,CGP,MODEL,KBASIS)
%
% Also pass in precomputed KBASIS structure.
%

%------------- start of initialization section -------------------------



% input arguments
fid = varargin{1};
cgp = varargin{2};
model = varargin{3};
Ncpus = cgp.Ncpus;

% if cgp.time
%     tic;
% end
if (nargin==3)
    kbasis = make_kbasis(fid,model,cgp.mat);
else
    kbasis = varargin{4};
end

% if (Ncpus > 1)
%     MPI_Data;    
%     preDist_temp.Ncpus = Ncpus;
%     if isstruct(kbasis)
%         p_kbasis = rmfield(kbasis,strvcat('ekxx','ekyy'));
%         p_kbasis.mat = 0;
%         preDist_temp.ix = kbasis.ix;
%         preDist_temp.iy = kbasis.iy;
%         for ii = 1:Ncpus
%             preDist_temp.ekxx = kbasis.ekxx(:,:,ii);
%             preDist_temp.ekyy = kbasis.ekyy(:,:,ii);
%             preDist{1} = preDist_temp;
%             MPI_Send(ii,'preDist');
%         end
%     else
%         p_kbasis.mat = 1;
%         for ii = 1:Ncpus
%             preDist_temp.ekxy = kbasis(:,:,ii);
%             preDist{1} = preDist_temp;
%             MPI_Send(ii,'preDist');
%         end
%         %clear kbasis;
%     end
%     p_kbasis.Ncpus = Ncpus;
%     clear preDist_temp preDist;
% end

% if isfield(model,'x') 
%     model = rmfield(model,strvcat('x','y','mask'));  % strip down model for efficient comm
% end

pv = model.pvec;

CC = zeros(1,200);  % set up variables for performance tracking
amp = zeros(model.reso,model.reso);
R2f = zeros(model.reso,model.reso,cgp.NIT);

% initialize the parameter array values as complex array
% pv((1:model.N)+model.N) = real(pv((1:model.N)+model.N))/model.wfr + 1i*imag(pv((1:model.N)+model.N))/model.wfi;

% bound frequency, and decay <=0
if isfield(cgp,'bounddecay')
    ind2 = 1+model.N:2*model.N;
    binz = false(1,model.N);
    pvBndDecay = abs(model.wfr*real(pv(ind2))) > cgp.bounddecay;
    pv([pvBndDecay binz]) = cgp.magdown*pv([pvBndDecay binz]);
    pv([binz pvBndDecay]) = 0;
    pv(ind2) = min(real(pv(ind2)),0) + 1i*imag(pv(ind2));
end

NF0 = 0;
niter = 1;
CCold = 1e+15;
CCnew = 1e+14;
modelcurr = model;

zeroed = 0;
nloop = 0;
if cgp.time, toc, fprintf('end of initialization\n'), end
%------------- end of initialization section ---------------------------   

if mod(niter,cgp.val) == 0,
    fprintf('%d %d %d %11.4e %11.4e\n',niter,cgp.NLIST(1),NF0,abs(CCold-CCnew),CCnew)
end
if mod(niter,cgp.time) == 0, toc, end

% begin iteration
tic;
for nlen = 1:cgp.NLENGTHS  % nested loops for progressive-length CG implementation
    fid.N = cgp.NLIST(nlen);  % current length
    fid.N = fid.N - mod(fid.N,Ncpus);
    modelcurr.pvec = pv;
    if (Ncpus < 2)
        fid.est = fidsyn(modelcurr,kbasis,fid.N);
    else
        fid.est = parallel_fidsyn(modelcurr,p_kbasis,fid.N);
    end
    %J0 = jeval(fid,modelcurr,cgp);

    while (abs(CCold-CCnew) > cgp.FLIST(nlen)*CCold || fid.N > NF0) && nloop < cgp.NIT
        nloop = nloop+1;  % keep track of iterations in each sub-loop
        
        if (Ncpus < 2)
            gradJ = s_grad(modelcurr,fid,kbasis,cgp);
        else
            gradJ = parallel_grad(modelcurr,fid,p_kbasis,cgp);
        end
        
        if mod(niter,cgp.time) == 0, toc, end
        % define direction of 1-D search for SD or CG (CG = 1);
        if fid.N ~= NF0
            g = -gradJ;
            gold = g;
            h = g;
        else
            gold = g;
            g = -gradJ;
            ga = max(0,real(((g-gold)*g')/(gold*gold')));
            h = g + ga*h;
        end
%         if mod(niter,cgp.vis) == 0
%             figure(13);subplot(2,1,1);title('gradients computed');
%             plot([real(h); imag(h)].');drawnow;
%         end

        modelcurr.pvec = pv;
        modelcurr.pvectmp = h;
        if (Ncpus < 2)
            delpf = Newton2(modelcurr,fid,kbasis,cgp);
        else
            delpf = parallel_Newton(modelcurr,fid,p_kbasis,cgp);
        end
        
        pv = pv + delpf*modelcurr.pvectmp;
        % bound frequency, and decay <=0
        if isfield(cgp,'bounddecay')
            pvBndDecay = abs(model.wfr*real(pv(ind2))) > cgp.bounddecay;
            pv([pvBndDecay binz]) = cgp.magdown*pv([pvBndDecay binz]);
            pv([binz pvBndDecay]) = 0;
            pv(ind2) = min(real(pv(ind2)),0) + 1i*imag(pv(ind2));
        end
        modelcurr.pvec = pv;
        if (Ncpus < 2)
            fid.est = fidsyn(modelcurr,kbasis,fid.N);
        else
            fid.est = parallel_fidsyn(modelcurr,p_kbasis,fid.N);
        end
        %fid.est = fidsyn(modelcurr,kbasis,fid.N);
        CC(niter) = jeval(fid,modelcurr,cgp);

        if mod(niter,cgp.plot) == 0
            figure(11); %subplot(2,1,2); 
            semilogy((CC(1:niter))); 
            title(strcat(num2str(fid.N),'  --  ',num2str(CC(niter)),'    ')); 
            drawnow
        end
        for nl = 1:model.N
            amp(model.ix(nl),model.iy(nl)) = pv(nl);
            expo1(model.ix(nl),model.iy(nl)) = pv(nl+model.N);
        end

        if (exist('fid.kmask'))
            % impose bandlimiting operator on result
            ampf(:,:,niter) = fftshift(ifft2(fftshift(fid.kmask.*fftshift(fft2(fftshift(amp))))));
        else
            ampf(:,:,niter) = amp;
        end
        expf(:,:,niter)=expo1;

        for nl = 1:model.N
            pv(nl) = ampf(model.ix(nl),model.iy(nl),niter);
        end

%         R2f(:,:,niter) = real(-model.wfr*squeeze(expf(:,:,niter)))/fid.delt;
        if mod(niter,cgp.vis) == 0
            figure(14);
            subplot(2,2,1); 
            imagesc(imag(squeeze(expf(:,:,niter))));axis image 
            title(['frequency ' num2str(fid.offr)]); ylabel('freq (Hz)');
            % subplot(2,2,2); 
            % plot(abs(squeeze(ampf(:,33,niter)))); 
            % title('abs(M0), col 33');
            subplot(2,2,3); 
            imagesc(-real(squeeze(expf(:,:,niter))));axis image 
            title('est. R2*'); ylabel('R2* (1/sec)');
            subplot(2,2,4); imagesc(abs(squeeze(ampf(:,:,niter)))); 
            axis image;colormap gray;
            subplot(2,2,4); 
            title(strcat('estimated magnitude, ',num2str(niter),' iterations')); 
            drawnow;
        end
        
        NF0 = fid.N;
        CCold = CCnew;
        CCnew = CC(niter);
        niter = niter + 1;
        if mod(niter,cgp.val) == 0, 
            fprintf('%d %d %d %11.4e %11.4e\n',niter,fid.N,NF0,abs(CCold-CCnew),CCnew)
        end
    end
end
%toc
model.pvec = pv;
model.CC = CC;
varargout{1} = model;
if nargout==2
    varargout{2} = kbasis;
end
%if cgp.time, toc, fprintf('end of CG algorithm\n'), end
